function g_ar = as_test_covariates(X, r0, con_dim)
% X: covariates to be conditioned on
% con_dim: the dimension of continuous variables that need to be transformed

n = size(X, 1);% number of observations
dx = size(X, 2);

if isempty(con_dim)
    X2 = X;
else
    % transform the continuous X
    Xc = X(:, con_dim);
    Xc2 = normcdf((Xc - mean(Xc))./std(Xc));    
%     sig_X = 1/n*(Xc - mean(Xc))'*(Xc - mean(Xc));
%     Xc2 = normcdf( (Xc - mean(Xc) )*chol(inv(sig_X + 0.0001*diag(diag(sig_X))))');

    X2 = X;
    X2(:, con_dim) = Xc2;
end

g_ar = cell(r0, 1);
% define the cube: cartesian product
for r = 1 : r0
    % divide each dimension to [0, 1/(2r)], [1/(2r), 2/(2r)], ...
    xr_grid = 0 : 1/(2*r) : 1;
    xr_grid(end) = xr_grid(end)+0.1; % to account for binary discrete variables

    % which part of the grid does an x fall into
    X2_grid_id = nan(n, dx);
    for ndx = 1 : dx
        X2_grid_id(:, ndx) = ...
            sum((X2(:, ndx)>=xr_grid(1:end-1) ...
            & X2(:, ndx)<xr_grid(2:end)).*(1 : 2*r),2);
    end

    % transform the grid assignments 
    [~, ~, ind_grid] = unique(X2_grid_id, 'rows');

    g_ar{r} = dummyvar(ind_grid);
end



